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Optimization  of  a Curve  Approximation 
Based  on  NURBS  interpolation 

Jerome  Lepine,  Frangois  Guibault, 
Marie-Gabrielle  Vallet,  and  Jean-Yves  Trepanier 


Abstract.  In  this  paper,  an  approach  is  presented  whereby  optimal 
spatial  positions  and  weights  of  a fixed  number  of  NURBS  control  points 
are  determined  using  a quasi-Newton  optimization  algorithm  in  order  to 
approximate  a general  planar  target  curve.  A method  for  constructing  an 
adequate  initial  solution  and  a valid  cost  function  based  on  interpolation 
error  are  introduced.  Convergence  of  the  iterative  process  is  assessed, 
and  the  final  interpolation  error  is  related  to  prescribed  manufacturing  or 
analysis  tolerances.  The  efficiency  of  the  approach  is  demonstrated  for 
actual  wing  profiles. 


§1.  Introduction 

The  problem  of  constructing  a cost  effective  approximation  of  a general  target 
curve  is  of  great  relevance  in  many  engineering  disciplines.  This  problem  has 
been  addressed  quite  thoroughly  in  the  context  of  polynomial  interpolation 
[1],  but  far  less  work  has  been  published  on  rational  approximation.  Indeed, 
weights  introduce  another  level  of  difficulty  in  the  theoretical  analysis  of  the 
approximation  error.  From  a practical  standpoint  though,  non-uniform  ra- 
tional B-Splines  [4]  (NURBS)  provide  more  degrees  of  freedom  for  a given 
number  of  control  points,  which  leads  naturally  to  smoother  curves. 

The  work  presented  here  introduces  a robust  numerical  approach  for  the 
determination  of  control  point  positions  and  weights  of  a NURBS  curve;  it 
can  be  used  to  construct  an  approximation  to  a general  target  planar  curve. 
In  the  context  of  wing  profile  design,  where  this  approach  has  been  applied 
[2,3],  very  significant  reductions  in  terms  of  data  size  and  noise  level  have  been 
observed. 

In  this  paper,  the  approximation  problem  is  first  presented,  and  the 
method  of  computation  of  the  approximation  error  discussed.  Next,  the  op- 
timization method  itself  is  presented,  including  the  choice  of  initial  solutions. 
Finally,  the  performance  of  the  method  is  evaluated  for  practical  test  cases, 
and  conclusions  are  drawn. 
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§2.  Approximation  Problem 
A NURBS  curve  is  defined  such  that 


A{u)  = YJRMPi  (1) 

t=0 


with 


Ri,P(u) 


NiiP{u)u>i 
o Nj,p(u)u)j 


(2) 


where  Pi  are  the  control  point  coordinates,  uj,  their  respective  weights,  Nl)P 
the  p-th  degree  B-spline  basis  functions  and  A{u)  the  position  of  a point  on 
the  curve.  The  basis  functions  are  obtained  through  a knot  vector,  which 
defines  the  functions’  break  points,  of  the  form 


p 


m — p-f  1 


Using  these  interpolation  functions,  the  problem  of  approximating  a general 
planar  curve  C(t)  can  be  stated  as  follows:  find  the  set  of  control  points  Pi 
and  weights  u>i  such  that  ||  A(u)  — C(t ) ||  is  minimized  in  a suitable  norm. 

Analytically,  the  L 2 norm  would  be  a natural  choice;  numerically,  though, 
for  a completely  general  taget  curve,  this  norm  can  only  be  approximated 
through  discretization.  Numerical  experiments  have  thus  been  carried  out  to 
develop  and  validate  a robust  computational  approach  for  the  determination 
of  the  appoximation  error.  Consideration  has  been  given  to  both  the  mean 
and  maximum  error,  as  well  as  to  the  level  of  continuity  of  the  target  curve. 
Three  classes  of  target  curves  have  been  considered:  curves  only  given  as  a 
set  of  points,  piecewise  linear  curves,  and  C 1 or  more  continuous  curves.  In 
all  cases,  the  mean  error 


Emea. 


k= 1 


(3) 


is  determined  by  summing  the  distance  (d*)  of  a set  of  points  chosen  on  the 
target  curve  to  their  respective  projection  on  the  approximation  curve,  and 
the  maximum  error 

Emax  — 77UZXi<fc<ndfc  (4) 

gives  the  largest  of  these  distances. 

As  a sample  of  these  experiments,  Fig.  1 illustrates  the  behaviour  of  the 
approximation  error  between  two  typical  target  curves,  a piecewise  linear  (a) 
and  a quadratically  interpolated  B-Spline  curve  (b),  and  their  approximation 
constructed  using  a NURBS  curve  with  13  control  points.  Both  target  curves 
were  specified  using  143  control  points.  As  can  be  observed  from  the  graphs  of 
Fig.  2,  a very  large  number  of  evenly  spaced  discretisation  points  must  be  used 
to  accurately  compute  both  the  mean  and  maximum  approximation  errors  for 
the  piecewise  linear  case.  The  behaviour  for  quadratic  test  cases  and  higher 
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Fig.  1.  Approximation  error  for  piecewise  linear  (left)  and  quadratic  (right) 
curve  (magnified  100X). 


Fig.  2.  Maximum  (left)  and  mean(right)  error  as  a function  of  the  number  of 
integration  points  for  the  piecewise  linear  target. 

degree  of  continuity  examples  (not  shown)  are  extremely  similar.  The  same 
graphs  of  Fig.  2 also  show,  as  a straight  line,  the  error  computed  using  only 
the  definition  control  points  of  the  target  curve. 

In  light  of  these  experiments,  it  was  determined  that  the  error  computed 
using  the  control  points  constitutes  an  adequate  bound  on  both  the  mean  and 
maximum  error  of  approximation,  and  it  can  be  computed  at  a fraction  of  the 
cost  of  using  evenly  spaced  discretization  points.  This  method  of  computing 
the  error  also  has  the  property  that  it  includes  naturally  the  case  of  target 
curves  given  as  a discrete  set  of  points,  which  is  not  a rare  case  in  many 
practical  applications. 


§3.  Optimization  Method 

Using  these  definitions  and  computational  method  of  the  approximation  er- 
ror, the  optimization  problem  can  be  further  specified  by  introducing  a cost 
function  of  the  form 

F(X)  = 2 X £mea  + E-max  i 

where  X is  the  vector  of  design  variables,  in  this  case  the  position  and  weights 
of  the  approximation  curve:  X = {xi,  y\,  wi,  X2,  ■ . . , xn,yn,  u>n}.  This  choice 
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of  a cost  function  significantly  accelerates  convergence  of  the  optimization 
process  by  including  both  the  maximum  error,  which  controls  the  quality  of 
the  final  approximation,  and  the  mean  error,  which  globally  compares  the 
quality  of  different  solutions. 

Clearly,  this  is  a non-linear  optimization  problem,  and  we  will  now  ex- 
amine the  chosen  solution  process,  including  the  choice  of  an  initial  solution. 

Solution  method 

The  primary  solution  method  used  was  the  second-order  quasi-Newton  me- 
thod, which,  given  a reasonably  close  initial  solution  Xo,  will  iteratively  con- 
verge towards  an  optimal  solution  using  the  relation 

Xk  + 1 = Xk  + Ofc  Sk , 

where  Sk  = —Hk  ■ VF(Xk)  is  the  direction  of  descent  vector,  and  a*,  the 
distance  of  descent  in  direction  Sk-  The  descent  vector  is  computed  using  the 
BFGS  [6]  algorithm,  based  on  a second  order  approximation  of  the  gradient 
of  F(X): 

VF(X)~VF(Xk)  + H(Xk)-6X, 

where  6X  = X — X k is  used  as  the  direction  of  descent  vector  ( Sk ).  Here 
H,  the  approximate  Hessian  matrix,  is  initially  set  to  identity  and  iteratively 
updated  using  the  relation 


Hk+i  = Hk  + 


Yk®Yk 

YkSk 


(Hk-Sk)®(Hk-Sk) 
Sk-  Hk-  Sk 


with  Yk  = VF(Xk+ i)  — VF(Xk)-  The  distance  of  descent  is  computed  using 
Armijo’s  rule  [5],  where  ak  = (5)™  and  m is  the  smallest  integer  such  that 
the  relationship 


F{Xk  + ak  Sk)  < F(Xk)  + aakVF(Xk)  • Sk 

with  a the  sufficient  descent  criterion,  which  must  be  chosen  between  0 and 
| (usually  set  to  10-4). 

Initial  solution 

In  most  cases,  the  optimization  method  described  above  will  find  a solu- 
tion, but  in  the  case  of  highly  non-linear  cost  functions  such  as  the  one  used  in 
this  problem,  it  is  impossible  to  determine  whether  the  minimum  found  is  the 
global  minimum  or  only  a local  one.  The  only  way  to  circumvent  this  difficulty 
is  to  proceed  with  many  optimizations,  and  select  the  best  minimum  as  the 
solution.  While  this  approach  could  be  unaffordable  if  no  clue  were  available 
about  the  solution,  it  can  be  implemented  relatively  cheaply  in  the  context  of 
curve  approximation,  where  many  good  initial  guesses  can  be  constructed. 

Specifically,  a set  of  initial  solutions  is  constructed  by  discretizing  the 
target  curve  using  a fixed  number  of  points  and  by  varying  the  concentration  of 
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points  along  the  curve.  Basically,  points  are  gathered  closer  together  in  regions 
of  high  curvature,  and  a shifting  constant  is  introduced  to  construct  various 
concentration  laws.  For  target  curves  of  continuity  less  than  C2,  curvature  is 
approximated  using  centered  differencing.  The  concentration  law  is  evaluated 
using 

Hu)  = ^jQ  C(v)+Ddv , 

where  C(v)  is  the  true  or  approximated  curvature  of  the  target  curve,  and  D 
the  shifting  constant. 

When  the  shifting  constant  becomes  large,  the  concentration  law  becomes 
almost  uniform.  In  practice,  sets  of  8 to  10  initial  solutions  are  constructed 
by  varying  D typically  between  1.0  and  10,  and  each  initial  solution  is  then 
optimized.  Fig.  3 shows  the  final  approximation  error  for  a run  where  D took 
the  values  {0.5, 1.0, 2.5, 3.0, 3.5, 4.5, 5.0, 6.0, 7.0}.  The  target  curve  for  this 
problem  is  a standard  NACA  2412  wing  profile,  and  9 control  points  are  used 
for  the  approximation,  which  leads  to  a 21  parameter  optimization  problem 
(the  two  endpoints  are  fixed).  Initial  weights  are  all  set  to  1.0. 

Fig.  3 vividly  illustrates  the  high  non-linearity  of  the  problem,  where  small 
variations  in  the  initial  solution  lead  to  completely  different  optimal  solutions, 
as  expressed,  for  example,  by  the  steep  variation  in  final  error  for  D = 4.5 
and  D --  5.0. 


§4.  Application 

We  will  now  look  at  how  this  approximation  method  performs  in  the  context 
of  a practical  application,  both  in  terms  of  data  reduction  and  approximation 
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Fig.  4.  Precision  of  the  approximation  as  the  number  of  control  points  increases. 

characteristics.  Two  aspects  of  the  approximation  are  of  particular  interest: 
precision  and  noise  level. 

The  application  consists  in  approximating  wing  profiles,  specified  either 
as  analytical  functions  or  experimental  sets  of  points.  In  this  context,  study 
[3]  of  the  combined  precision  levels  dictated  by  both  manufacturing  tolerances 
and  precision  for  analysis  purposes  indicates  that  a precision  of  the  order  of 
8 X 10~5  is  sufficient. 

Precision 

Fig.  4 illustrates  the  evolution  of  emax  as  the  number  of  control  points  of 
the  approximation  curve  is  increased.  Again  the  target  curve  is  the  NACA 
2412  profile.  As  can  be  observed,  the  increase  in  precision  of  the  approximation 
is  very  regular  when  8 control  points  or  more  are  used.  In  this  case,  the 
required  precision  of  8 x 10-5  is  obtained  with  only  9 control  points.  Extensive 
experiments  [3]  involving  numerous  types  of  wing  profiles  have  shown  that  the 
required  level  of  precision  can  almost  always  be  attained  with  13  control  points 
or  less.  These  numbers  have  to  be  compared  with  the  number  of  points  needed 
to  discretely  represent  a profile  with  the  same  precision,  which  can  be  shown 
to  be  of  the  order  of  150.  The  approximation  method  thus  offers  excellent 
control  over  the  precision  of  the  resulting  curve,  while  reducing  by  more  than 
an  order  of  magnitude  the  amount  of  data  used  for  representation. 

Experiments  have  also  been  carried  out  in  order  to  determine  whether 
the  introduction  of  weights  in  the  formulation  had  an  impact  on  precision. 
Similar  precision  tests  carried  out  using  B-Splines  instead  of  NURBS  have 
shown  that  exactly  1.5  more  control  points  were  required  for  B-Splines  to 
obtain  comparable  levels  of  precision.  This  increase  in  the  number  of  control 
points,  however,  has  a significant  impact  on  noise  level. 

Noise 

For  many  engineering  applications,  such  as  airplane  wing  design,  noise 
level  is  often  a bigger  concern  than  absolute  precision  level.  In  that  respect, 
the  NURBS  approximation  method  performs  remarkably  well,  mainly  because 
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Fig.  5.  Curvature  of  the  Boeing  A8  profile  and  of  its  approximation  near  the 
leading  edge. 

of  the  small  number  of  control  points  needed.  Noise  appears  as  small  fluctua- 
tions in  the  curvature  of  a curve,  particularly  when  control  points  are  gathered 
closely  together  in  regions  of  high  curvature.  As  shown  in  Fig.  5,  the  proposed 
approximation  method  can  significantly  reduce  noise  in  cases  where  the  target 
curve  presents  important  fluctuations.  Of  course,  this  reduction  of  the  noise 
level  can  only  be  accomplished  as  a trade-off  to  the  precision  of  the  approx- 
imation. For  example,  the  precision  of  the  13  control  point  approximation 
of  the  Boeing  A8  profile  of  Fig.  5 is  9.2  x 10~5,  which  is  slightly  above  the 
usual  tolerance  level  for  this  application;  better  precision  could  be  obtained 
by  including  a few  more  control  points,  but  this  would  inevitably  introduce 
more  noise. 


§5.  Conclusion 

We  have  presented  a method  of  approximation  for  a general  planar  curve  that 
permits  a significant  reduction  in  the  size  of  data  and  garantees  a desired  level 
of  precision.  The  main  advantages  of  this  approach  are 

• generality, 

• full  automation, 

• low  noise. 

By  varying  a single  parameter,  the  shifting  constant  in  the  construction  of 
the  concentration  law,  as  many  initial  solutions  as  needed  are  generated  for 
a given  number  of  points.  Each  solution  is  optimized  independently,  and  the 
solution  with  minimum  error  is  kept.  Because  of  the  typically  small  number 
of  control  points  required  - of  the  order  of  10  to  15  - very  smooth  curves 
are  obtained,  which  is  a very  important  characteristic  for  many  engineering 
applications. 

Because  of  the  significant  reduction  in  the  number  of  free  parameters 
used  to  represent  a curve,  the  approximation  method  is  now  being  used  as 
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a first  step  in  a shape  optimization  procedure  of  wing  profiles.  This  work  is 
also  currently  being  extended  to  three  dimensional  cases,  where  the  method 
is  now  used  to  approximate  wing  surfaces. 
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